function VARbs = AE_doCholSVARbootstrap_single(VAR,nboot,clevel)
% Code based on Gertler and Karadi (2015)
% Wild Bootstrap: uses original estimates and draws (with replacement) the residuals.


jj=1;
% Removes a constant
res = detrend(VAR.res,'constant');

% Bootstrap simulations (indexed by jj)
while jj<nboot+1

  % T x 1 dimension array of -1 or 1 random numbers with prob 0.5-0.5 (Rademacher distribution)
  rr = 1-2*(rand(VAR.T,1)>0.5);

  % T x n randomly choose the sign of the time T shocks (all)
  resb = (res.*(rr*ones(1,VAR.n)))';

  % Preallocate matrix of endogenous variables
  varsb = zeros(VAR.p+VAR.T,VAR.n);

  % Initial values
  varsb(1:VAR.p,:)=VAR.vars(1:VAR.p,:);

  % Use baseline estimates to generate endogenous variables, using newly simulated shocks
  for j=VAR.p+1:VAR.p+VAR.T
    lvars      = (varsb(j-1:-1:j-VAR.p,:))'; % vector of lagged variables
    varsb(j,:) = lvars(:)'*VAR.bet(1:VAR.p*VAR.n,:) + VAR.bet(VAR.p*VAR.n+1,:) + resb(:,j-VAR.p)';
  end

  VARBS      = VAR;
  VARBS.vars = varsb;

  % Estimate again the VAR, and Cholesky (computes IRF)
  VARBS = AD_doCholSVAR_single(VARBS);

  % Save results into VARbs (output of function)
  for j=1:VAR.k;
    irs = VARBS.irs(:,:,j);
    VARbs.irs(:,jj,j) = irs(:);
  end
                 
  jj=jj+1;   

end

% Compute confidence intervals
for j=1:VAR.k  
  VARbs.irsL(:,:,j)=reshape(quantile(VARbs.irs(:,:,j)',(1-clevel/100)/2),VAR.irhor,size(VARBS.irs,2));
  VARbs.irsH(:,:,j)=reshape(quantile(VARbs.irs(:,:,j)',1-(1-clevel/100)/2),VAR.irhor,size(VARBS.irs,2));
end


